###
library(clubSandwich)
library(lmtest)
library(sandwich)
### Large Cities
Pult3 <- readRDS("Pult3.rds")
summary(Pult3$Civilian)
BigOnes <- lm(Employees ~ LRP + FreeSpace  + Rugged +
                `Avg Height (ft)` + I(`Avg Height (ft)`^2) +
                `Land Area (sq mi)` + I(`Land Area (sq mi)`^2) + Pop+
                factor(Year) +
                factor(`NBC Affiliate (City, Call Sign)`),
              data = Pult3[which(Pult3$Year>Pult3$YearStarted),])
###
cluster_var <- Pult3$`NBC Affiliate (City, Call Sign)`[Pult3$Year > Pult3$YearStarted]

cl_vcov <- vcovCL(BigOnes, cluster = cluster_var)

# Coefficients with clustered SEs
coefs <- coeftest(BigOnes, vcov = cl_vcov)
coefs_filtered <- coefs[!grepl("^factor\\(Year\\)|^factor\\(`NBC Affiliate", rownames(coefs)), , drop = FALSE]
coefs_filtered
BigEmp <- coefs_filtered["LRP",]
# Print the result
print(coefs_filtered)
####
BigOnes2 <- lm(Civilian ~ LRP + FreeSpace  + Rugged +
                 `Avg Height (ft)` + I(`Avg Height (ft)`^2) +
                 `Land Area (sq mi)` + I(`Land Area (sq mi)`^2) + Pop+
                 factor(Year) +
                 factor(`NBC Affiliate (City, Call Sign)`),
               data = Pult3[which(Pult3$Year>Pult3$YearStarted),])
###

cluster_var <- Pult3$`NBC Affiliate (City, Call Sign)`[Pult3$Year > Pult3$YearStarted]

cl_vcov <- vcovCL(BigOnes2, cluster = cluster_var)

# Coefficients with clustered SEs
coefs <- coeftest(BigOnes2, vcov = cl_vcov)
coefs_filtered <- coefs[!grepl("^factor\\(Year\\)|^factor\\(`NBC Affiliate", rownames(coefs)), , drop = FALSE]
CivEmp <- coefs_filtered["LRP",]
# Print the result
print(coefs_filtered)
###
Pult3$Treated <-0
Pult3$Treated[which(Pult3$Year>Pult3$YearStarted)] <- 1
DnDL<- lm(Employed ~ Treated+ Pop+
            factor(Year) +
            factor(ID),
          data = Pult3)

cluster_var <- Pult3$`NBC Affiliate (City, Call Sign)`

# Cluster-robust variance-covariance matrix
cl_vcov <- vcovCL(DnDL, cluster = cluster_var)

# Print clustered SEs
Large <- coeftest(DnDL, vcov = cl_vcov)
LargeTWFE <- Large["Treated",]
Large[c(2,3),]
DnDLC <-  lm(Civilian ~ Treated+
               factor(Year) + Pop+
               factor(ID),
             data = Pult3)
cluster_var <- Pult3$`NBC Affiliate (City, Call Sign)`

# Cluster-robust variance-covariance matrix
cl_vcov <- vcovCL(DnDLC, cluster = cluster_var)

# Print clustered SEs
CivT <- coeftest(DnDLC, vcov = cl_vcov)
CivT[c(2,9),]
CivTWFE <- CivT["Treated",]
###
#saveRDS(results,"results.rds")

results <- readRDS("results.rds")
# load in this data MERGED with expenditure data

MainExp <- lm(PoliceExp ~ LRP + FreeSpace  + Rugged +
                `Avg Height (ft)` + I(`Avg Height (ft)`^2) +
                Area + I(Area^2)+
                factor(Year.x) +
                factor(`NBC Affiliate (City, Call Sign)`)+Pop+Expend,
              data = results[results$Year.x > results$YearStarted,])
vcov_clustered <- vcovCL(MainExp, cluster = results$`NBC Affiliate (City, Call Sign)`[results$Year.x > results$YearStarted])

# Coefficient test with clustered SEs
coef_table <- coeftest(MainExp, vcov = vcov_clustered)

# Remove factor terms
coef_table_clean <- coef_table[!grepl("^factor\\(", rownames(coef_table)), ]

# Print the cleaned table
print(coef_table_clean)
PolExp <- coef_table_clean["LRP",]
###

MainExp <- lm(EducExpend ~ LRP + FreeSpace  + Rugged +
                `Avg Height (ft)` + I(`Avg Height (ft)`^2) +
                Area + I(Area^2)+
                factor(Year.x) +
                factor(`NBC Affiliate (City, Call Sign)`)+Pop+Expend,
              data = results[results$Year.x > results$YearStarted,])
vcov_clustered <- vcovCL(MainExp, cluster = results$`NBC Affiliate (City, Call Sign)`[results$Year.x > results$YearStarted])

# Coefficient test with clustered SEs
coef_table <- coeftest(MainExp, vcov = vcov_clustered)

# Remove factor terms
coef_table_clean <- coef_table[!grepl("^factor\\(", rownames(coef_table)), ]

# Print the cleaned table
print(coef_table_clean)
EducExp <- coef_table["LRP",]
####

MainExp <- lm(Expend ~ LRP + FreeSpace  + Rugged +
                `Avg Height (ft)` + I(`Avg Height (ft)`^2) +
                Area + I(Area^2)+
                factor(Year.x) +
                factor(`NBC Affiliate (City, Call Sign)`)+Pop+EducExpend,
              data = results[results$Year.x > results$YearStarted,])
vcov_clustered <- vcovCL(MainExp, cluster = results$`NBC Affiliate (City, Call Sign)`[results$Year.x > results$YearStarted])

# Coefficient test with clustered SEs
coef_table <- coeftest(MainExp, vcov = vcov_clustered)

# Remove factor terms
coef_table_clean <- coef_table[!grepl("^factor\\(", rownames(coef_table)), ]

# Print the cleaned table
print(coef_table_clean)
HwyExp <- coef_table["LRP",]
###
### All Cities
#saveRDS(stacked_data2,"stacked_data2.rds")
stacked_data2 <- readRDS("stacked_data2.rds")

stacked_data2$Rugged2 <- stacked_data2$Max_Height_ft - stacked_data2$Min_Height_ft
stacked_data2$NBCAffiliateCityCallSign <- recode(
  stacked_data2$NBCAffiliateCityCallSign,
  "San Francisco, KNTV" = "KRON",
  "KCRA-TV (Sacramento)" = "KCRA",
  "KSBW (Monterey)" = "KSBW",
  "KSBW (Salinas)" = "KSBW",
  "Monterey, KSBW" = "KSBW",
  "KRGV-TV" = "KRGV",
  "KFDX-TV" = "KFDX",
  "WOAI-TV" = "WOAI",
  "WBRC/WABT" = "WBRC",
  "wdtn" = "WDTN",
  "wdtn / cleve" = "WDTN",
  "wkyc" = "WKYC",
  "wlio" = "WLIO",
  "wlwt" = "WLWT",
  "wfmj" = "WFMJ",
  "WGAl" = "WGAL",
  "WVVA (Bluefield)" = "WVVA"
)
MainGrad <- lm(Emp ~ LRP + FreeSpace + Rugged +
                 Avg_Height_ft + I(Avg_Height_ft^2) +
                 Land_Area + I(Land_Area^2) + Population +
                 factor(Year) +
                 factor(NBCAffiliateCityCallSign),
               data = stacked_data2[which(stacked_data2$Year > stacked_data2$year_start & stacked_data2$Year > 1950),])

ctest <-coeftest(MainGrad, vcov = vcovCL(MainGrad, cluster = stacked_data2$NBCAffiliateCityCallSign[which(stacked_data2$Year > stacked_data2$year_start & stacked_data2$Year > 1950)]))

###
vars_to_keep <- c("LRP", "FreeSpace", "Rugged", "Avg_Height_ft", "I(Avg_Height_ft^2)", 
                  "Land_Area", "I(Land_Area^2)", "Population","(Intercept)")

ctest[rownames(ctest) %in% vars_to_keep, , drop = FALSE]
AllEmp <- ctest["LRP",]
####
### Graphic
expenditure_df <- tibble::tibble(
  Outcome = c("Education Expenditure", "Highway Expenditure", "Police Expenditure"),
  Estimate = c(EducExp[1], HwyExp[1], PolExp[1]),
  SE = c(EducExp[2], HwyExp[2], PolExp[2])
) %>%
  mutate(
    CI_lower = Estimate - 1.96 * SE,
    CI_upper = Estimate + 1.96 * SE
  )

# Plot
max_abs <- max(abs(c(expenditure_df$CI_lower, expenditure_df$CI_upper)))
ggplot(expenditure_df, aes(x = Estimate, y = Outcome)) +
  geom_point(size = 3) +
  geom_errorbarh(aes(xmin = CI_lower, xmax = CI_upper), height = 0.2) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_x_continuous(limits = c(-max_abs, max_abs)) +
  labs(
    x = "Expenditure (in $1000s)",
    y = "Outcome Variable",
    title = "Effect of Signal Strength Loss on Expenditure"
  ) +
  theme_minimal(base_size = 14)
### Employment Graph

data <- data.frame(
  Category = c("Total Officers Employed (All Cities)", 
               "Total Officers Employed (Large Cities)", 
               "Civilian Police Employees (Large Cities)"),
  Coefficient = c(AllEmp[1], BigEmp[1], CivEmp[1]),
  StdError = c(AllEmp[2], BigEmp[2], CivEmp[2])
)
data$CI_Lower <- data$Coefficient - 1.96 * data$StdError
data$CI_Upper <- data$Coefficient + 1.96 * data$StdError
data$Category <- factor(data$Category, levels = c(
  "Total Officers Employed (All Cities)", 
  "Total Officers Employed (Large Cities)", 
  "Civilian Police Employees (Large Cities)"
))

# Calculate confidence intervals (95%)
data$CI_Lower <- data$Coefficient - 1.96 * data$StdError
data$CI_Upper <- data$Coefficient + 1.96 * data$StdError

# Find y-axis limits to center zero
y_limit <- max(abs(c(data$CI_Lower, data$CI_Upper))) + 0.05  # Adding a bit of padding

# Plot
data$Category <- factor(data$Category, levels = c(
  "Civilian Police Employees (Large Cities)", 
  "Total Officers Employed (Large Cities)", 
  "Total Officers Employed (All Cities)"
))

# Calculate confidence intervals (95%)
data$CI_Lower <- data$Coefficient - 1.96 * data$StdError
data$CI_Upper <- data$Coefficient + 1.96 * data$StdError

# Find y-axis limits to center zero
y_limit <- max(abs(c(data$CI_Lower, data$CI_Upper))) + 0.05  # Adding a bit of padding

# Plot
ggplot(data, aes(x = Category, y = Coefficient)) +
  geom_point(size = 3, color = "black") +  # Black points for estimates
  geom_errorbar(aes(ymin = CI_Lower, ymax = CI_Upper), width = 0.2, color = "black") +  # Black error bars
  geom_hline(yintercept = 0, linetype = "dotted", color = "black", linewidth = 1) +  # Dotted zero line
  theme_minimal() +
  labs(title = "Regression Coefficients for Reduced Form IV with 95% CIs",
       x = "Outcome",
       y = "Coefficient Regressing LRP on Law Enforcement Employees") +
  coord_flip() +  
  theme(
    legend.position = "none",  # Remove legend
    text = element_text(color = "black"),  # Ensures black text
    plot.title = element_text(hjust = 0.5)  # Centers title
  ) +
  ylim(-y_limit, y_limit)
###
### Two Way Fixed Effects
# Need to load in data WITH 1950
# This was done post the others because 50 is in a weird format
#saveRDS(stacked_data3,"stacked_data3.rds")
stacked_data3 <- readRDS("stacked_data3.rds")


table(stacked_data3$Year)


stacked_data3$NBCAffiliateCityCallSign <- recode(
  stacked_data3$NBCAffiliateCityCallSign,
  "San Francisco, KNTV" = "KRON",
  "KCRA-TV (Sacramento)" = "KCRA",
  "KSBW (Monterey)" = "KSBW",
  "KSBW (Salinas)" = "KSBW",
  "Monterey, KSBW" = "KSBW",
  "KRGV-TV" = "KRGV",
  "KFDX-TV" = "KFDX",
  "WOAI-TV" = "WOAI",
  "WBRC/WABT" = "WBRC",
  "wdtn" = "WDTN",
  "wdtn / cleve" = "WDTN",
  "wkyc" = "WKYC",
  "wlio" = "WLIO",
  "wlwt" = "WLWT",
  "wfmj" = "WFMJ",
  "WGAl" = "WGAL",
  "WVVA (Bluefield)" = "WVVA"
)



DnD <- lm(Emp~Treated+factor(Year) +
            factor(ID)+Population,data=stacked_data3)
#DnD <- lm(Emp~Treated+factor(Year) +
#            factor(ID)+LPop,data=stacked_data3)

clustered_se <- vcovCL(DnD, cluster = stacked_data3$NBCAffiliateCityCallSign)
full_results <- coeftest(DnD, vcov = clustered_se)
filtered_results <- full_results[rownames(full_results) %in% c("Treated", "Population"), ]
filtered_results
AllTWFE <- filtered_results["Treated",]
###

##
results$Treated <- 0
results$Treated[which(results$Year.x>results$YearStarted)] <- 1

DnD <- lm(PoliceExp~Treated+factor(Year.x) +
            factor(ID)+Pop+EducExpend,data=results)

clustered_se <- vcovCL(DnD, cluster = results$`NBC Affiliate (City, Call Sign)`)
full_results <- coeftest(DnD, vcov = clustered_se)
filtered_results <- full_results[rownames(full_results) %in% c("Treated", "Pop"), ]
filtered_results
ExpTWFE <- filtered_results["Treated",]
###
# Create the dataset
data_employment <- data.frame(
  Category = c("Police Employees (All Cities)", "Officers Employed (Large Cities)", "Civilians Employed (Large Cities)"),
  Coefficient = c(AllTWFE[1], LargeTWFE[1], CivTWFE[1]),
  StdError = c(AllTWFE[2], LargeTWFE[2], CivTWFE[2])
)

# Set default z-score to 1.96 (for 95%)
z_scores <- rep(1.96, nrow(data_employment))

#


# Calculate confidence intervals
data_employment$CI_Lower <- data_employment$Coefficient - z_scores * data_employment$StdError
data_employment$CI_Upper <- data_employment$Coefficient + z_scores * data_employment$StdError
# Employment Plot
ggplot(data_employment, aes(x = Category, y = Coefficient)) +
  geom_point(size = 3, color = "black") +  # Points for estimates
  geom_errorbar(aes(ymin = CI_Lower, ymax = CI_Upper), width = 0.2, color = "black") +  # Confidence intervals
  geom_hline(yintercept = 0, linetype = "dotted", color = "black", linewidth = 1) +  # Dotted zero line
  theme_minimal() +
  labs(title = "Effect of Dragnet Introduction on Police Employment",
       x = "Outcome Variable",
       y = "Change in Number of Employees") +
  coord_flip() +
  theme(legend.position = "none")
### Expenditure Plot
    
exp_df <- tibble::tibble(
  Outcome = "Police Expenditures",
  Estimate = ExpTWFE[1],
  SE = ExpTWFE[2] 
) %>%
  mutate(
    CI_lower = Estimate - 1.96 * SE,
    CI_upper = Estimate + 1.96 * SE
  )

# Plot
ggplot(exp_df, aes(x = Estimate, y = Outcome)) +
  geom_point(size = 3) +
  geom_errorbarh(aes(xmin = CI_lower, xmax = CI_upper), height = 0.2) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  labs(
    x = "Effect in Dollars per 1000 People",
    y = "Outcome Variable",
    title = "Effect of Dragnet Introduction on Police Expenditure"
  ) +
  xlim(-100, 150) +  # Match original plot scale
  theme_minimal(base_size = 14)

###


